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The affect of demographic stochasticity of a system of globally coupled chaotic maps is considered. 
A two-step model is studied, where the intra-patch chaotic dynamics is followed by a migration step 
that coupled all patches; the equilibrium number of agents on each site, N, controls the strength 
of the discreteness-induced fluctuations. For small N (large fluctuations) a period-doubling cascade 
appears as the coupling (migration) increases. As N grows an extremely slow dynamic emerges, 
leading to a flow along a one-dimensional family of almost period 2 solutions. This manifold be- 
come a true solutions in the deterministic limit. The degeneracy between different attractors that 
characterizes the clustering phase of the deterministic system is thus the N — > oo limit of the slow 
dynamics manifold. 
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The dynamics of coupled chaotic maps have attracted 
a lot of interest in the last decades, following the pio- 
neering works of Kaneko [TJ [2] . A substantial part of the 
study is focused around the paradigmatic model of glob- 
ally coupled maps, where many fundamental results like 
mutual synchronization, dynamical clustering and glassy 
behavior were demonstrated [5] . The universal character 
of the chaotic dynamics makes the coupled maps model 
relevant to many phenomena, ranging from neural sys- 
tems and human body rhythms to coupled lasers and 
cryptography [3J. 

Here we consider the effect of demographic stochastic- 
ity (shot noise) on various phases of a globally coupled 
system. This problem emerges naturally while applying 
the theory to spatially extended ecologies. 

Many old [5] and recent [B] experiments suggest that 
the well-mixed dynamics of simple ecosystems (single 
species or victim-exploiter system) are extinction-prone, 
and that the system acquires stability only due to its 
spatial structure, a result supported also by numerical 
simulations of many models pj |5]. The extended (spa- 
tial) system survives due to the possibility of migration 
among spatial patches. This migration should be large 
enough to allow for recolonization of empty patches be 
emigrants. On the other hand [9], too much migration is 
also dangerous, as it leads to global synchronization, in 
which case the system acts essentially as a single, well- 
mixed patch, with its vulnerability to extinction. 

The globally coupled system, which obeys, 

si +1 = (i-y)F(si)+^F( s i), a) 

is a natural and popular framework to discuss the dynam- 
ics of so-called meta-populations [10 on various patches 
with migration between the patches [5]. Here, Sj is the 
population density on the i-th site and F is the chaotic 
map. v is the migration parameter (the chance of an in- 
dividual agent to leave its site) and 1 < i < L, where L is 
the number of patches. To address the problem of extinc- 
tion properly, however, it is necessary to account for the 



discrete nature of the population and the absorbing char- 
acter of the zero population state. This is achieved by 
the addition of demographic (shot) noise to the system. 
Clearly, if the local populations arc all large, this effect is 
tiny. However, for small and moderate populations, the 
effects as we shall see can be very important and give rise 
to new phenomena. This is somewhat surprising, since 
it is widely assumed that chaos generates its own noise, 
and the addition of other noise should not induce a qual- 
itative change in the dynamics. However, one important 
feature of coupled maps is the appearance of attractive 
regular orbits for certain value of the coupling; it turns 
out that once the system is kicked from these orbits it 
follows a long excursion on its way back (this issue will 
be discussed in a separate publication). The presence 
of continuous noise may thus change drastically the ob- 
served dynamics. The results of the deterministic theory 
(like clustering) should reappear, though, in the large N 
limit, where N is the number of particles at a patch. 

Let us begin our discussion by specifying a stochastic 
model. Again we have a collection of L sites, where the 
local population density s\ is now replaced by an integer 
n\. The update proceeds in two steps. First, the re- 
production and competition generates a new value of n, . 
This value is taken to be drawn from a Poisson distribu- 
tion with mean F(rii), where F(n) is a chaotic map. In 
this paper, we take as our choice of F the paradigmatic 
Packer map [11], F(n) — rne~ n / N , where r is the growth 
factor. Also, we have fixed the value of r = 20, which is 
well in the chaotic regime of the deterministic map. (We 
have chosen the Ricker dynamics just because it sim- 
plifies our numerics; our results hold for a broad range 
of different maps, and we believe that any chaotic map, 
including the logistic one, should exhibit similar behav- 
ior.) The second phase is the dispersal phase, in which 
with probability v, each of the inhabitants of every site 
can decide to leave and pick a new site at random. This 
model is essentially similar to that used by Hamilton and 
May [12] to study optimal dispersal rates, except for the 
chaotic nature of the on-site reproduction/competition 
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dynamics of the present model. 

We have simulated this system directly using Monte- 
Carlo technique for L — 10000. However for the purpose 
of analysis, it is more convenient to study, as do Hamilton 
and May, the L — > oo limit. This limit is completely char- 
acterized by a probability distribution ip^, the chance of 
a given site to have n individuals at time t. The probabil- 
ity of having m individuals after birth and competition is 
fixed by n, and the probability distribution for the num- 
ber of individuals leaving to other sites is then fixed by 
to. Since we are dealing with an infinite reservoir, the 
probability distribution for the number of incoming in- 
dividuals is fixed by A, the average population after the 
birth/competition phase, which in turn is fixed by ip'- 



X = F{m)= J2^F(m). 



m=0 



Accordingly, t/j dynamics follows, 
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where fJ,(m) — F(m)(\ — v) + v\. The update rule is a 
linear transformation of the probability vector ip which 
conserves probability. The transformation matrix A4 = 
M(X) is a functional of the input state, since it depends 
on A, which depends of It is this dependence on the 
input state that renders the problem nonlinear, and gives 
it its rich structure. 

We consider first a constant (period 1) solution. Pick- 
ing an arbitrary value for A one gets a Markov matrix 
A4(X) from Due to the Markov property this matrix 
must admit an invariant eigenvector (a right eigenvector 
with eigenvalue 1) ip n (X). This procedure is consistent 
iff V-'n(^) an d A satisfy Eq. This consistency con- 

dition determines A for the period 1 solution, which in 
turn determines ?/> of the period 1 solution. A period 2 
solution is obtained in a similar way: Picking arbitrary 
values for Ai and A 2 the Markov matrix M(\2)M(\i) 
must admit an invariant eigenvector i\) n . The solution is 
consistent iff the system satisfies the two auxiliary con- 
ditions Ax = 5>n-F(n) and A 2 = £ [M(\i)ip] n F(n). 
Extending this procedure one may find orbits of higher 
periodicity by searching through the space of quartets 
and octets of A-s with the appropriate auxiliary condi- 
tions. 

One trivial solution always exists for period 1 orbits: 
the absorbing state for which A = and tp n — <5„o- It 
turns out that for small values of v this is the only solu- 
tion and the probability vector converges exponentially 
quickly to S n0 (see Figure [T] where the results for strong 
stochasticity, N = 5, are summarized). This is as ex- 
pected, since there is a finite probability of an individual 
site to go extinct and without sufficient dispersal to en- 
able recolonization, more and more sites go extinct as 
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FIG. 1: Upper: Solution branches for N = 5, showing the 
period 1, 2 and 4 branches as a function of v. Lower: Prob- 
ability Distributions ip n for the period 1 solution at v — 0.2 
and for the two phases of the period 2 solution at v = 0.25. 
The points are connected to guide the eye. 



time goes on. Above a critical value of v the system 
quickly adopts a nontrivial period- 1 configuration. The 
function n{y) is plotted in Figure [l] The corresponding 
distribution ij) n is characterized by two peaks, as exem- 
plified in the lower panel of [l] It shows two peaks, each 
of which has half the probability. The system decom- 
poses into two clusters that oscillates 180° out of phase 
with respect to the other. Since the two clusters have 
equal weight, the overall occupancy of the system is time- 
independent. 

What about the stability of the period- 1 solution? The 
above numerical technique of course only identifies solu- 
tions and does not say anything about their stability. 
One can prove, however, that any periodic orbit must 
become unstable as v — > 1. As pointed out by Durrett 
and Levin [7], in that case the occupation of a site just 
before the reaction step is a Poisson distribution with a 
mean given by the total population in the last step. As 
a result, n satisfies the iterative map, 
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where F is the deterministic map. In the Ricker case 
the resulting map for n is also unimodal and the result- 
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ing transformation is chaotic in the regime of parameters 
considered here. Thus any periodic orbit must lose its 
stability as v approaches unity. 

All the considerations and the results presented so far 
are general and are independent of the strength of the 
stochasticity, which is inversely proportional to N. We 
turn now to consider the differences between strong and 
weak stochasticity and the semi-deterministic (large N) 
limit. 

Let us refer again to Figure [T] where the results for 
N = 5 (strong stochasticity) are graphed. There is 
a range of v for which the period 1 is stable, as can 
be verified by MC simulation. Increasing v this time- 
independent state goes unstable and the system under- 
goes a forward (supercritical) bifurcation to a period-two 
state. This state losses stability in favor of a period-4 so- 
lution and so on; our results suggest that in that case 
a cascade of period doubling supercritical bifurcations 
emerges until, at some v, there is a transition to global 
chaos. The period two solution also supports a bimodal 
distribution, but now the weight of the two peaks are 
not equal so the overall population takes different value 
as the peaks switch positions. A similar scenario happens 
in the period 4 regime. 

Comparing our results with the deterministic system 
discussed in [3], in the presence of strong demographic 
noise one may observe the fully synchronized phase (when 
v approaches one and the map shows global chaos) but 
what about the turbulent phase? This phase appears in 
the deterministic system when the coupling is too small 
and different patches oscillate incoherently. Clearly the 
extinction phase is a consequence of turbulence. How- 
ever, the transition from the turbulent phase to the clus- 
ter phase in the deterministic system occurs at much 
higher v than the extinction transition for the stochastic 
system (which is exponentially small in N). The answer 
to this puzzle lies in the fact that there is no clear distinc- 
tion between the cluster phase and the turbulent phase 
in the stochastic system. As N or v are reduced, the dis- 
tribution broadens and the identifiable peaks wash out. 
Furthermore the total signal in both phases is time inde- 
pendent. Only in the infinite N limit can a sharp distinc- 
tion between the two phases be made. In this limit, the 
distribution function becomes two delta function peaks 
at infinite N in the cluster phase, and have finite support 
in the turbulent phase.. 

Let us now focus our attention on the clustering phase 
of the deterministic system. . One of the features of this 
phase is observed here: the system clusters spontaneously 
into two peaks and the global population follows a pe- 
riodic orbit. However, there are two important differ- 
ences: first, in the deterministic system for the same v 
many possible solutions exist, each corresponds to dif- 
ferent height ratio between the peaks, while here for any 
migration parameter only one stable solution survives, so 
the infinite degeneracy that characterizes the determinis- 



tic dynamics disappears. Second, the sharp (delta) peaks 
of the deterministic solutions are replaced by smooth dis- 
tributions. One expects, however, that the stochastic 
system converges to the deterministic one as N increases. 
How does that happen? 

Let us start to increase N. For N — 10 (results not 
shown) the bifurcation from period 1 to period 2 is back- 
ward (subcritical) , while the bifurcation from period 2 to 
period 4 is still forward even at N = 20. The location 
of the bifurcation from period 1 to 2 is almost indepen- 
dent of N, but the bifurcation from 2 to 4 is strongly N 
dependent, moving to smaller v as N increases. In fact, 
by N — 40 it has already moved to the backward branch 
of the period 2 solution. This situation is summarized in 
the upper panel of Fig. 2, where the period 1, 2 and 4 
solution branches are traced out for TV = 60. 

The most important issue is how the deterministic con- 
tinuum of period two solutions is recovered as N grows 
to infinity. As explained above the period 2 solutions are 
identified by searching for all pairs of Ai, A2 that admit an 
invariant eigenvector ip such that Ai = ^ipnF^n), and 
A 2 = J2 F ( n ) [• / ^(^i)V']„- It turns out that, for large N, 
there is a range of Ai, A2 for which the Markov matrix 
M21 = -M(A 2 ).M(Ai) admits, in addition to its invariant 
eigenvector, an additional eigenstate ip with an eigen- 
value very close to 1, say, 1 — e. Thus, up to a small term 
(for N — 80, e.g., e = 10 -12 ) any linear combination of 
the first and the second eigenvectors imitates the real in- 
variant state until t ~ 1/e. Within this time horizon one 
has, effectively, a continuous family of invariant eigenvec- 
tors of M21, (1 — a)tp + aip- The two auxiliary conditions 
no longer are sufficient to determine a solution. We call 
these solutions for which we ignore e a quasi- solution, of 
which there exists a continuous family depending on a. It 
turns out that e decreases sharply with increasing N; the 
deterministic limit emerges from this continuous family 
of solutions as explained below. 

In Fig. [3] we see all this exemplified in a simulation, 
where we plotted Tit as a function of t, for N = 60, 
v = 0.21. We see that for times less than roughly 5 • 10 5 , 
the system exhibits an essentially period 2 type behavior, 
with an extremely slow drift of the two states. Suddenly, 
beyond this point, the system converts to a period 4 be- 
havior. A good way of analyzing the drift is to plot ri2t+i 
vs. n,2t, as seen in Fig. [2] (lower panel). If the system 
had a true period 2 orbit, this graph would show a single 
point. Instead, the drift converts this into a curve. The 
points on this curve coincide precisely with the above de- 
scribed quasi-solutions, a number of which are indicated 
by circles. The true solutions are represented by triangles 
(period 2) and a diamond (period 1) in the figure. The 
system drifts to larger amplitude oscillations, until the 
instability is encountered and the it goes to a period 4 
orbit, represented by two dots in the figure. While there 
are quasi-solutions (as well as a true solution a = 0) be- 
yond this point, they are not dynamically relevant due 
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FIG. 2: Upper: Solution branches for N = 60, showing the 
period 1, 2 and 4 branches as a function of v. Lower: The 
return map nzt+i versus n,2t, indicated by black dots, for N = 
60, v — 0.21, taken from the simulation depicted in Fig. [3] 
The trajectory - a drift away from the green triangle, then 
split to the period 4 orbit - is determined by the line of quasi 
solutions indicated by the red circles. The triangles represent 
true period 2 solutions for which a — 0, and the period 1 
solution is indicated by a red diamond. The inset is a blowup 
of the relevant section of the upper panel, indicating the true 
solutions and the slow flow through the quasi-solutions. The 
solutions for v = 0.21 are marked by the same symbols as 
in the main figure. For these parameters the dynamics close 
to the green triangle is so slow that measuring the dynamics 
becomes impractical. When TV — ■> oo, all the region between 
the green triangle and the splitting point (indicated by full 
green line in the inset) becomes marginally stable. 

to their strong instability. 

In summary, then, we have seen that adding demo- 
graphic noise to a globally coupled chaotic map has a 
marked effect on the dynamics, leading in fact to very 
regular dynamics for intermediate coupling strength. As 
the noise strength is reduced, there appears an exponen- 
tially long scale, which goes over to the continuous family 
of solutions seen in the no noise limit. 

Beyond the problems considered here, our work sug- 
gests a novel mechanism for stochastic-deterministic 
("quantum-classical") correspondence. In general it is 
assumed that the infinite number of solutions that char- 
acterizes the deterministic case emerge from the discrete 
eigenstates of the stochastic theory when the level spac- 



400 
300 


1 


1 

J 


200 






100 











1 


i 



2xl0 5 4xl0 5 6xl0 5 

t 



FIG. 3: The average population n t versus t for N = 60, 
v = 0.21. Data obtained by direct integration of the Mas- 
ter equation. 



ing approaches zero in the weak stochasticity limit. Here 
we observed a different scenario, where only two eigen- 
states provide us with a continuum of deterministic so- 
lutions at the N — > oo limit. 
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